Python notebook¶
Met deze Notebook genereren we de volgende output:
vertical_resistance_layer1: een nieuwe verticale weerstand voor de deklaagwinter\peil_laag#_#.idf: nieuwe peilen voor de winter in de verschillende (sub) lagenzomer\peil_laag#_#.idf: nieuwe peilen voor de zomer in de verschillende (sub) lagen
Een aantal relevante tussenproducten:
peilen.gpkg: de geaggregeerde vaste/zomer/winter en flexibele peilen naar 1 zomer en 1 winterpeil, wanneer deze kan worden berekendwinterpeil_ruw.tif,winterpeil.tif,winterpeil.idf: het ruwe en geinterpoleerde resultaat van de verrastering van de kolomwinterpeilinpeilen.gpkgzomerpeil_ruw.tif,zomerpeil.tif,zomerpeil.idf: het ruwe en geinterpoleerde resultaat van de verrastering van de kolomzomerpeilinpeilen.gpkg
In deze Notebook worden alle tussen en eindproducten stap-voor-stap gegenereerd met duiding en Python-code.
try: # to bypass ImportError at rasterio.__init__ from rasterio.__version ....
import rasterio
except ImportError:
from osgeo import gdal
import rasterio
from dask import array # If you don't do this imod.idf.open() (version 0.9.0) will give an AttributeError: module 'dask' has no attribute 'array'
import geopandas as gpd
import imod
import numpy as np
from pathlib import Path
import shutil
from rasterio.features import rasterize
from rasterio.fill import fillnodata
import xarray as xr
from nhi_data import Ondergrond
import warnings
warnings.filterwarnings("ignore") # To avoid logging RuntimeWarning in rasterio.fill.fillnodata
Specificaties iMODFLOW inlezen¶
We lezen de specificatie van het iMODFLOW model in. We gebruiken daarvoor de laag PEIL_LAAG1_1.IDF, waarvan we de ruimtelijke attributen overnemen.
oppervlaktewater_dir = Path(r"d:\projecten\D2304.HDSR_TEO_debieten\02.model\data\hydromedah\oppervlaktewater")
template_idf = oppervlaktewater_dir.joinpath("winter", "PEIL_LAAG1_1.IDF")
template_da = imod.idf.open(template_idf)
shape = template_da.shape
height, width = template_da.shape
transform = imod.util.transform(template_da)
_, xmin, xmax, _, ymin, ymax = imod.util.spatial_reference(template_da)
Vervangen deklaagweerstand¶
Hier vervangen we de verticale weerstand van de deklaag uit het LHM4.1
Inladen NHI ondergrond¶
We maken een Ondergrond-object, voor de bounding-box van het HDSR MODFLOW model. We zoeken de juiste laag op de NHI-geoserver
ondergrond = Ondergrond(bbox=(xmin,ymin,xmax,ymax))
layer = ondergrond.get_layers(filter="resistance_layer1")[0]
print(layer)
ondergrond2:vertical_resistance_layer1_lhm41
Download en conversie naar GeoTIFF en IDF¶
Met dit object kunnen we de laag downloaden en opslaan als GeoTiff en IDF
ondergrond.download_layer(layer, "vertical_resistance_layer1.tif")
ondergrond.download_layer(layer, "vertical_resistance_layer1.idf")
Vervangen peil en peilgebieden¶
Hier vervangen we de oppervlaktewaterpeilen van het iMODFLOW model
Inladen peilgebieden¶
De peilgebieden zijn op 04-07-2023 gedownload (gegenereerd op 23 jun 2023 11:04) via de HDSR Infovijver als GeoJSON.
peilengebieden_gdf = gpd.read_file("Peilgebieden.geojson", engine="pyogrio")
peilengebieden_gdf.explore()
Aggregeren peilen naar zomer- en winterpeil¶
We lezen de geometrieen in van de peilgebieden en kopieren het attribuut SOORTSTREEFPEIL (soort_streefpeil) met de volgende enumeratie:
1: vast peil2: flexibel peil3: zomer- en winterpeil4,5,99: geen peil vastgesteld/peil onbekend
We willen de volgende kolommen vullen:
code: waterschapscodenodata: geen zomer- winterpeil bekend =True, andersFalsebetrouwbaar_peil:winterpeil<=zomerpeil=True, andersFalsezomerpeilwinterpeil
We projecteren de geometriën naar Rijksdriehoekstelsel (epsg 28992) en vullen de kolommen soort_streefpeil, code uit de oorspronkelijke data en zetten de waarden van de kolom nodata op False.
peilen_gdf = gpd.GeoDataFrame(peilengebieden_gdf["geometry"])
peilen_gdf.to_crs(28992, inplace=True)
peilen_gdf["soort_streefpeil"] = peilengebieden_gdf["SOORTSTREEFPEIL"]
peilen_gdf["code"] = peilengebieden_gdf["CODE"]
peilen_gdf["nodata"] = False
Vast Peil¶
Wanneer VASTPEIL == -999 of niet ingevuld, dan geven we dit aan in nodata. Voor de peilgebieden waarvan gegevens bekend zijn geldt zomerpeil = winterpeil = vastpeil
peil_mask = peilengebieden_gdf["SOORTSTREEFPEIL"] == 1
# VASTPEIL = -999 komt voor, dit maskeren we en leggen we vast als no data
data_mask = (peilengebieden_gdf["VASTPEIL"] != -999) & ~peilengebieden_gdf["VASTPEIL"].isna()
peilen_gdf.loc[(peil_mask & ~data_mask), ["nodata"]] = True
# vullen zomerpeil en winterpeil
peilen_gdf.loc[(peil_mask & data_mask), ["zomerpeil"]] = peilengebieden_gdf[(peil_mask & data_mask)]["VASTPEIL"]
peilen_gdf.loc[(peil_mask & data_mask), ["winterpeil"]] = peilengebieden_gdf[(peil_mask & data_mask)]["VASTPEIL"]
Flexibel peil¶
In zowel ONDERPEIL als BOVENPEIL komt -999 voor, maar niet gelijktijdig. In deze gevallen stellen we ONDERPEIL gelijk aan BOVENPEIL of vise-versa. Vervolgens brekenen we zomerpeil en winterpeil:
zomerpeil = winterpeil = (ONDERPEIL + BOVENPEIL) / 2
peil_mask = peilengebieden_gdf["SOORTSTREEFPEIL"] == 2
# bij -999 in bovenpeil stellen we bovenpeil gelijk aan onderpeil
nodata_mask = peil_mask & (peilengebieden_gdf["BOVENPEIL"] == -999)
peilengebieden_gdf.loc[nodata_mask, ["BOVENPEIL"]] = peilengebieden_gdf[nodata_mask]["ONDERPEIL"]
# bij -999 in onderpeil stellen we onderpeil gelijk aan bovenpeil
nodata_mask = peil_mask & (peilengebieden_gdf["ONDERPEIL"] == -999)
peilengebieden_gdf.loc[nodata_mask, ["ONDERPEIL"]] = peilengebieden_gdf[nodata_mask]["BOVENPEIL"]
peilen_gdf.loc[peil_mask, ["zomerpeil"]] = (peilengebieden_gdf[peil_mask]["ONDERPEIL"] + peilengebieden_gdf[peil_mask]["BOVENPEIL"]) / 2
peilen_gdf.loc[peil_mask, ["winterpeil"]] = peilen_gdf[peil_mask]["zomerpeil"]
Zomer- en winterpeil¶
Hier komen geen ontbrekende, danwel -999 waarden voor, dus deze nemen we direct over
peil_mask = peilengebieden_gdf["SOORTSTREEFPEIL"] == 3
peilen_gdf.loc[peil_mask, ["zomerpeil"]] = peilengebieden_gdf[peil_mask]["ZOMERPEIL"]
peilen_gdf.loc[peil_mask, ["winterpeil"]] = peilengebieden_gdf[peil_mask]["WINTERPEIL"]
Onbrekende peilen¶
Wanneer peilen ontbreken, geven we dat aan door nodata op True te zetten
peil_mask = peilengebieden_gdf["SOORTSTREEFPEIL"] >= 4
peilen_gdf.loc[peil_mask, ["nodata"]] = True
Indicatie betrouwbaarheid¶
Voor alle gebieden met zomerpeil en winterpeil controleren we of winterpeil kleiner of gelijk aan zomerpeil is.
data_mask = ~peilen_gdf["nodata"]
peilen_gdf.loc[data_mask, ["peil_betrouwbaar"]] = peilen_gdf[data_mask]["winterpeil"] <= peilen_gdf[data_mask]["zomerpeil"]
Wegschrijven¶
peilen_gdf.to_file("peilen.gpkg")
Conversie naar imod¶
Het resultaat in peilen_gdf wordt weggeschreven naar idf-bestanden. Van de verrasterde peilen schrijven we tevens GeoTifs weg, zodat deze in GIS kunnen worden geladen.
crs = peilen_gdf.crs
profile={ # profiel voor GeoTifs
"driver": "GTiff",
"dtype":rasterio.dtypes.float32,
"nodata": -999,
"width" :width,
"height": height,
"count": 1,
"crs": crs,
"transform": transform}
Data mask¶
We zorgen dat de Utrechtse Heuvelrug (PG2179) en Nederrijn en Lek (PG2112) uit de peilen-rasters blijven. We houden daar de oorspronkelijke waarden. Het deel waar we peilen voor verrasteren schrijven we weg als mask.tif.
shapes = (
(geom, 1)
for geom in peilen_gdf[~peilen_gdf["code"].isin(["PG2179", "PG2112"])]["geometry"]
)
mask_array = rasterize(
shapes=shapes,
out_shape=shape,
transform=transform,
all_touched=True)
with rasterio.open("mask.tif", "w+", **profile) as dst:
dst.write(mask_array, 1)
Verwerking in IDFs¶
Voor de seizoenen zomer en winter voeren we de volgende stappen uit:
- Verrasteren van de zomer en winterpeilen waar deze data bevatten. We schrijven deze weg in
zomerpeil_ruw.tifenwinterpeil_ruw.tif - Het opvullen van
nodatacellen binnen het mask door middel van interpolatie. We schrijven deze weg alszomerpeil.tifenwinterpeil.tif - Het overschrijven van de waarden in alle
PEIL_LAAG*.IDFbestanden uit het oorsprokelijke model met de regels:nodatain de oorspronkelijke peilen blijftnodatain de nieuwe peilen- Wanneer er een er een bodemhoogte is wordt het nieuwe peil gelijk of hoger dan de bodemhoogte
for seizoen in ["zomer", "winter"]:
# 1. Verrasteren van de zomer en winterpeilen waar deze data bevatten. We schrijven deze weg in `zomerpeil_origineel.tif` en `winterpeil_origineel.tif`
shapes = (
(geom, value)
for geom, value in zip(peilen_gdf[~peilen_gdf.nodata]["geometry"], peilen_gdf[~peilen_gdf.nodata][f"{seizoen}peil"])
)
peil_array = rasterize(
shapes=shapes,
out_shape=shape,
fill=-999,
transform=transform,
all_touched=False)
with rasterio.open(f"{seizoen}peil_ruw.tif", "w+", **profile) as dst:
dst.write(peil_array, 1)
# 2. Het opvullen van `nodata` cellen binnen het mask door middel van interpolatie. We schrijven deze weg als `zomerpeil.tif` en `winterpeil.tif
mask = np.where((peil_array == -999) & (mask_array == 1), 0, 1)
mask = np.where(mask_array == 0, np.NaN, mask)
peil_array = fillnodata(
peil_array,
mask=mask,
max_search_distance=20,
smoothing_iterations=0
)
peil_array = np.where(mask_array == 1, peil_array, -999)
with rasterio.open(f"{seizoen}peil.tif", "w+", **profile) as dst:
dst.write(peil_array, 1)
dataset = imod.rasterio.open(f"{seizoen}peil.tif")
imod.idf.save(f"{seizoen}peil", dataset)
# 3. Het overschrijven van de waarden in alle `PEIL_LAAG*.IDF` bestanden
out_dir = Path(seizoen)
if out_dir.exists():
shutil.rmtree(out_dir)
out_dir.mkdir()
datasets = imod.idf.open_dataset(oppervlaktewater_dir.joinpath(seizoen,"PEIL_LAAG*.IDF"))
for laag, peil_da in datasets.items():
# - `nodata` in de oorspronkelijke peilen blijft `nodata`in de nieuwe peilen
peil_da = xr.where(
np.isnan(peil_da),
peil_da,
np.where(peil_array == -999, peil_da, peil_array)
)
# - Wanneer er een er een bodemhoogte is wordt het nieuwe peil gelijk of hoger dan de bodemhoogte
bodemhoogte_name = f"BODEMHOOGTE_LAAG{laag[len(laag)-3:]}"
bodemhoogte_da = imod.idf.open(
oppervlaktewater_dir.joinpath(seizoen,f"{bodemhoogte_name}.IDF")
)
peil_da = xr.where(
np.isnan(peil_da),
peil_da,
xr.where(
np.isnan(bodemhoogte_da),
peil_da,
xr.where(
bodemhoogte_da < peil_da,
peil_da,
bodemhoogte_da
) # zowel bodemhoogte als peil zijn niet NaN
)
)
# Wegschrijven IDF
imod.idf.write(out_dir / f"{laag}.idf", peil_da)